/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunriset is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Implementation of methods for the thermal_equilibrium_grain class.

// $Id$

#include "config.h"
#include "thermal_equilibrium_grain.h"
#include "interpolatort.h"
#include "threadlocal.h"
#include "vecops.h"
#include <cassert>

#ifdef WITH_CUDA
#include "cuda_grain_temp.h"
#endif

using namespace mcrx;
using namespace blitz;
using namespace std;

mcrx::thermal_equilibrium_grain::thermal_equilibrium_grain
(const std::string& file_name, const Preferences& p,
 T_float minsize, T_float maxsize) :
  emitting_grain(file_name, minsize, maxsize), s_(*this),
  accuracy_(1e-3)
{
  n_threads_ = p.getValue("n_threads", int(1), true);
  bind_threads_ = p.getValue("bind_threads", false, true);
  bank_size_ = p.getValue("thread_bank_size", int(n_threads_), true);

  accuracy_ = p.getValue("grain_temp_accuracy", T_float(1e-3), true);
  use_lookup_ = p.getValue("use_grain_temp_lookup", bool(false), true);
#ifdef WITH_CUDA
  use_cuda_ = p.getValue("use_cuda", bool(false), true);
#else
  use_cuda_ = false;
#endif

  // now call our own setup statically
  thermal_equilibrium_grain::setup();
}

void mcrx::thermal_equilibrium_grain::setup()
{
  // call base class setup first. this is safe even if we call from
  // the derived constructor and it's already been called, because it
  // can be called several times.
  grain_data::setup();
  temp_table_updated_=false;
}


void mcrx::thermal_equilibrium_grain::setup_interpolators() const
{    
  // set up temp interpolators
  temp_interpolators_.resize(n_size());

  cout << "Setting up temperature interpolation tables" << endl;
  const int npoints=200;
  const T_float low_temp=3.0;
  const T_float hi_temp=1500;
      s_.setverbose(true);

  // loop over sizes
  for (int s = 0;s<sizes().size(); ++s) {
    // First find limits of the interpolation table corresponding to
    // temp limits for this grain size. We get this from -func when
    // absorption_=0
    current_grain_= s;
    heating_ = 0;
    const T_float low_abs = log10(-func(low_temp));
    assert(low_abs==low_abs);
    const T_float hi_abs = log10(-func(hi_temp));
    assert(hi_abs==hi_abs);
    vector<T_float> absorption_points;
    for(int j=0; j<npoints; ++j) {
      const T_float abs=pow(10., (hi_abs-low_abs)*j/(npoints-1)+low_abs);
      absorption_points.push_back(abs);
    }
    temp_interpolators_[s].setAxis(0, absorption_points);

    array_1 temps(absorption_points.size());
    // loop over temp points
    for (int j=0; j<absorption_points.size(); ++j)  {
      heating_ = absorption_points[j];
      T_float xx = log10(heating_) - 
	1.39*log10(asizes()(s)) + 
	0.126*pow(log10(asizes()(s)),2);
      T_float guess =
	pow(10,1.86+0.189*xx+3.41e-3*xx*xx);
  
      guess = (guess<5)?5:guess;
      temps(j) = solve_for_T(s, heating_, guess);
    }
    ASSERT_ALL(temps==temps);
    temp_interpolators_[s].setPoints(temps);
  }
  temp_table_updated_=true;
}


/** Solves for the grain temperature for size bin s, starting with the
    specified guess. The absorption vector should be loaded into
    absorption_ before this function is called. */
mcrx::T_float
mcrx::thermal_equilibrium_grain::solve_for_T(int s, T_float h,
					     T_float guess) const
{
  // the temp calculation has to be done in SI
  assert(units().get("wavelength")=="m");
  assert(units().get("length")=="m");
  assert(units().get("luminosity")=="W");

  if(h==0) 
    return 0;

  current_grain_ = s;
  heating_ = h;
  // set the accuracy for energy conservation
  s_.setaccuracy(accuracy_*h, accuracy_);
  int c=0;
  const T_solver::returntype solution = s_.solve(guess);
  const T_float T0 = solution.value;
  if(solution.status) {
    cerr << "Fatal: Failed to find temp solution for absorption " << h << " @T=" << T0 << endl;
    throw 1;
  } 
   
  DEBUG(3,cout << "TEMP FUNCTION " << s << '\t' << heating_ << '\t'	\
	<< T0 << endl;);

  assert(T0>=0);
  return T0;
}


T_float 
mcrx::thermal_equilibrium_grain::
interpolate_T(int s, T_float heating) const
{
  assert(heating>=0);
  assert(heating<1e300);

  // make sure temp interpolators are updated
  if(!temp_table_updated_)
    setup_interpolators();

  T_float T;
  const T_float lower = temp_interpolators_[s].lower_range()[0];
  const T_float upper = temp_interpolators_[s].upper_range()[0];
  if (heating >lower ) {
    assert(heating > lower);
    if(heating < upper) {
      assert(heating < upper);
      T = temp_interpolators_[s].interpolate(heating);
    }
    else {
      // abs>upper range, solve with upper temp range as guess
      assert(heating > upper);
      try {
	T = solve_for_T(s, heating, temp_interpolators_[s].interpolate(upper));
      }
      catch(...) {
	cerr << "Failed solving for temperature for size " << s << ". Interpolator range is " << lower << " - " << upper << endl;
	throw;
      }
    }
  }
  else {
    // T<lower range. Let's in that case just set it to 0K and be
    // done with it. There is so little emission anyway and the peak
    // is outside our wavelength range.
    assert(heating < lower);
    T = 0;
  }

  assert(T>=0);

  return T;
}
